// Bibliotecas
#include <stdio.h>
#include <cstdio>
#include <iostream>
#include <cstdlib>
#include <string>
#include <fstream>
#include <sstream>
#include <vector>
#include <cmath>
#include <sys/time.h>

//25
#define Nthreads 25

using namespace std;

__global__ void div_per_diag(float* a, float* b, int J_ORDER) {
    float elem;
    // Calculate the column index of the Pd element, denote by x
    int x; //= threadIdx.x + blockIdx.x * blockDim.x;
    // Calculate the row index of the Pd element, denote by y
    int y; //= threadIdx.y + blockIdx.y * blockDim.y;

    //int seq = x + y * gridDim.x * blockDim.x;
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;

    y = seq / J_ORDER;
    x = seq - J_ORDER * y;

    if ((x < J_ORDER) && (y < J_ORDER)) {
        if (x != y) {
            elem = a[y * J_ORDER + y];
            a[y * J_ORDER + x] /= elem;
            b[y] /= elem;
        }
    }
}

__global__ void zero_diag(float* a, int J_ORDER) {
// Calculate the column index of the Pd element, denote by x
    int x; //= threadIdx.x + blockIdx.x * blockDim.x;
    // Calculate the row index of the Pd element, denote by y
    int y; //= threadIdx.y + blockIdx.y * blockDim.y;

    //int seq = x + y * gridDim.x * blockDim.x;
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;

    y = seq / J_ORDER;
    x = seq - J_ORDER * y;

    if ((x < J_ORDER) && (y < J_ORDER)) {
        if (x == y) {
            a[y * J_ORDER + y] = 0;
        }
    }
}

__global__ void sum_row(float *a, int *test, int J_ORDER) {
    // Calculate the column index of the Pd element, denote by x

    //int seq = x + y * gridDim.x * blockDim.x;
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;

    if ((seq < J_ORDER)) {
        float sum = 0.0f;
        for (int i = 0; i < J_ORDER; i++)
            sum += a[seq * J_ORDER + i];
        if (sum > 1.0f) {
            *test = 1;
        }
    }

}

__global__ void d_copy(float *dx, float *b, int J_ORDER) {
    // Calculate the column index of the Pd element, denote by x
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    if ((seq < J_ORDER)) {
        dx[seq] = b[seq];
    }
}

__global__ void d_solve(float *s, float *dx, float *b, float *a, float *d,
        int J_ORDER) {
    // Calculate the column index of the Pd element, denote by x
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    if ((seq < J_ORDER)) {
        float sum = .0f;
        for (int j = 0; j < J_ORDER; j++) {
            if (j != seq)
                sum += a[seq * J_ORDER + j] * dx[j]; // pode ser invertido
        }
        //printf("x %d sum %f\n", x, -sum);
        s[seq] = b[seq] - sum;
        d[seq] = abs(s[seq] - dx[seq]);
        //if(x == J_ORDER-1) printf(">> %f\n", s[x]);
    }
}
__global__ void d_abs(float *dx, int J_ORDER) {
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    if ((seq < J_ORDER)) {

        dx[seq] = abs(dx[seq]);
    }

}

__global__ void d_max(float *dx, float *max, int J_ORDER) {
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    float out = .0f;
    if (seq == 0) {
        for (int i = 0; i < J_ORDER; i++) {
            out = out < dx[i] ? dx[i] : out;
        }
    }
    *max = out;
}

__global__ void d_ver(float *a, int i, float *dx, float *o, int J_ORDER) {
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    if ((seq < J_ORDER)) {
        dx[seq] = a[i * J_ORDER + seq] * dx[seq];
    }

}
__global__ void d_ver1(float *a, int i, float *dx, float *o, int J_ORDER) {
    int seq = (threadIdx.x + threadIdx.y * blockDim.x)
            + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.x;
    if (0 == seq) {
        float sum = .0f;
        for (int j = 0; j < J_ORDER; j++) {
            sum += dx[j];
        }
        *o = sum;
    }
}

bool test_converg(float *dev_A, int J_ORDER) {
    int *dev_test;
    int test = 0;

    cudaMalloc((void**) &dev_test, sizeof(int));
    cudaMemcpy(dev_test, &test, sizeof(int), cudaMemcpyHostToDevice);

    // kernel invocation code
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    sum_row<<<dimGrid, dimBlock>>>(dev_A, dev_test, J_ORDER);

    cudaMemcpy(&test, dev_test, sizeof(int), cudaMemcpyDeviceToHost);

    cudaFree(dev_test);
    return test > 0 ? false : true;
}

void copy(float *x, float *b, int J_ORDER) {
    // vector para interacoes
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    d_copy<<<dimGrid, dimBlock>>>(x, b, J_ORDER);
}

void solve(float *s, float *x, float *b, float *a, float *d, int J_ORDER) {
    // vector para interacoes
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    d_solve<<<dimGrid, dimBlock>>>(s, x, b, a, d, J_ORDER);

}

float max(float *x, int J_ORDER) {
    float max = .0f;
    float *dev_max;

    cudaMalloc((void**) &dev_max, sizeof(float));
    cudaMemcpy(dev_max, &max, sizeof(float), cudaMemcpyHostToDevice);

    // kernel invocation code
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    d_max<<<dimGrid, dimBlock>>>(x, dev_max, J_ORDER);

    cudaMemcpy(&max, dev_max, sizeof(float), cudaMemcpyDeviceToHost);

    cudaFree(dev_max);
    return max;
}

float erro(float *x, float *d, int J_ORDER) {
    // vector para interacoes
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    d_abs<<<dimGrid, dimBlock>>>(x, J_ORDER);

    return (float) (max(d, J_ORDER) / max(x, J_ORDER));
}

float ver(float *a, int J_ROW_TEST, float *x, int J_ORDER) {
    float value;
    float *dev_value;

    cudaMalloc((void**) &dev_value, sizeof(float));

    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    d_ver<<<dimGrid, dimBlock>>>(a, J_ROW_TEST, x, dev_value, J_ORDER);
    d_ver1<<<dimGrid, dimBlock>>>(a, J_ROW_TEST, x, dev_value, J_ORDER);
    cudaMemcpy(&value, dev_value, sizeof(float), cudaMemcpyDeviceToHost);

    cudaFree(dev_value);

    return value;
}

void jacobi_richardison(float *A, float *B, int J_ORDER, int J_ROW_TEST,
        float J_ERROR, int J_ITE_MAX) {
    int size1 = J_ORDER * sizeof(float);
    int size2 = J_ORDER * J_ORDER * sizeof(float);

    // variaveis host
    float *NA = (float*) malloc(size2);
    float *NB = (float*) malloc(size1);

    // variaveis no dispositivo
    float *dev_A;
    float *dev_B;
    float *dev_D;
    float *dev_X;
    float *dev_S;

    // allocate memory on the GPU
    cudaMalloc((void**) &dev_A, size2);
    cudaMalloc((void**) &dev_B, size1);

    // transfer M and N to device memory
    cudaMemcpy(dev_A, A, size2, cudaMemcpyHostToDevice);
    cudaMemcpy(dev_B, B, size1, cudaMemcpyHostToDevice);

    // kernel invocation code
    dim3 dimBlock(Nthreads, Nthreads);
    dim3 dimGrid(J_ORDER / Nthreads, J_ORDER / Nthreads);
    div_per_diag<<<dimGrid, dimBlock>>>(dev_A, dev_B, J_ORDER);

    zero_diag<<<dimGrid, dimBlock>>>(dev_A, J_ORDER);

    // verifica convergencia
    if (!test_converg(dev_A, J_ORDER)) {
        cudaFree(dev_A);
        cudaFree(dev_B);
        //printf("passou\n");
        return;
    }

    // aloca a memoria na GPU
    cudaMalloc((void**) &dev_X, size1);
    cudaMalloc((void**) &dev_D, size1);
    cudaMalloc((void**) &dev_S, size1);

    int k;
    copy(dev_X, dev_B, J_ORDER);
    for (k = 0; k < J_ITE_MAX; k++) {
        solve(dev_S, dev_X, dev_B, dev_A, dev_D, J_ORDER);
        copy(dev_X, dev_S, J_ORDER);
        if (erro(dev_S, dev_D, J_ORDER) <= J_ERROR)
            break;
    }

    cudaMemcpy(dev_A, A, size2, cudaMemcpyHostToDevice);

    cout << "Iterations: " << k << endl;
    cout << "ROW_TEST: " << J_ROW_TEST << " => ";
    cout << "[" << ver(dev_A, J_ROW_TEST, dev_X, J_ORDER) << "] aprox "
            << B[J_ROW_TEST] << endl;

    // transfer dev_C e dev_B from device
    cudaMemcpy(NA, dev_A, size2, cudaMemcpyDeviceToHost);
    cudaMemcpy(NB, dev_X, size1, cudaMemcpyDeviceToHost);
    /*
     cout << "matriz A* e matriz *B" << endl;
     for(int i = 0; i < 1; i++) {
     for(int j = 0; j < J_ORDER; j++) {
     cout << NA[i*J_ORDER+j] << " ";
     }
     cout << endl;
     }

     for(int j = 0; j < J_ORDER; j++) {
     cout << NB[j] << " ";
     }
     cout << endl;
     */

    free(NA);
    free(NB);
    cudaFree(dev_A);
    cudaFree(dev_B);
    cudaFree(dev_D);
    cudaFree(dev_S);
    cudaFree(dev_X);
}

int main(int argc, char* argv[]) {
    static int J_ORDER;
    static int J_ROW_TEST;
    static float J_ERROR;
    static int J_ITE_MAX;

    // interadores
    int i, j;

    // variavel para pegar as linhas do arquivo que esta em fin
    string line;
    ifstream fin(argv[1]); // ifstream de arquivo

    // Leitura do arquivo
    if (fin.is_open()) {
        string buffer;
        getline(fin, line, '\n'); // pega a linha 1 --> J_ORDER
        buffer = line;

        getline(fin, line, '\n'); // pega a linha 2 --> J_ROW_TEST
        buffer += " " + line;

        getline(fin, line, '\n'); // pega a linha 3 --> J_ERROR
        buffer += " " + line;

        getline(fin, line, '\n'); // pega a linha 4 --> J_ITE_MAX
        buffer += " " + line;

        istringstream iss(buffer);

        //DESCARREGANDO O BUFFER
        iss >> J_ORDER;
        iss >> J_ROW_TEST;
        iss >> J_ERROR;
        iss >> J_ITE_MAX;

        // Matrizes do sistema linear
        float *A, *B;
        A = (float*) malloc(J_ORDER * J_ORDER * sizeof(float));
        B = (float*) malloc(J_ORDER * sizeof(float));

        // Matriz A
        for (i = 0; i < J_ORDER; i++) {
            getline(fin, line, '\n');
            istringstream iss_line(line);

            for (j = 0; j < J_ORDER; j++) {
                iss_line >> A[i * J_ORDER + j];
            }
        }

        // Matriz B
        buffer.clear();
        for (j = 0; j < J_ORDER; j++) {
            getline(fin, line, '\n');
            buffer += " " + line;
        }

        istringstream iss_line(buffer);
        for (j = 0; j < J_ORDER; j++) {
            iss_line >> B[j];
        }

        // Chamada do metodo numerico
        jacobi_richardison(A, B, J_ORDER, J_ROW_TEST, J_ERROR, J_ITE_MAX);

        free(A);
        free(B);

        // fecha arquivo
        fin.close();

    } else
        cout << "not found " << argv[1] << endl; // mensagem de erro do arquivo

    return 0;
}

